/*
    IceCode toolkit - maths
    Copyright (C) 2011 Alfredo 'IceCoder' Mungo

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef _ICECODE_PYTHAGORAS_MATRIX
#define _ICECODE_PYTHAGORAS_MATRIX

#include <initializer_list>
#include "Exception.hpp"

namespace icecode {
  namespace maths {

    /*@
      DESCRIPTION:
       A generic matrix class.

      ARGUMENTS:
       M:
        Number of rows.
       N:
        Number of columns.
       T:
        Element type.
    */
    template<int M, int N = M, typename T = double>
    class Matrix
    {
      static_assert(M > 0, "The matrix must have positive number of rows.");
      static_assert(N > 0, "The matrix must have positive number of columns.");
      typedef Matrix<M, N, T> MTYPE;

    public:

      /*
	Creates an identity matrix.
       */
      Matrix() { identity(); }

      Matrix(std::initializer_list<T> args) {
	if(args.size() == N * M)
	  {
	    typename std::initializer_list<T>::iterator it = args.begin();
	    int i = 0;
	    do{
	      data[i++] = *it;
	      it++;
	    } while(it != args.end());
	  }
	else throw PYTHAGORAS_INVALID_DIMENSION_EXCEPTION;
      }

      Matrix(const Matrix& m) { *this = m; }

      bool operator ==(MTYPE& m)
      {
	for(int i = 0; i < N * M; i++)
	  if(data[i] != m.data[i])
	    return false;

	return true;
      }

      bool operator != (MTYPE& m)
      {
	return !(*this == m);
      }

      /*
	Returns a row.
      */
      T* operator [](int index)
      {
	return &data[index * M];
      }

      template <int L> //L=columns of 'm' matrix
      Matrix<M, L, T> operator *(Matrix<N, L, T>& m)
      {
	Matrix<M, L, T> res;

	for(int i = 0; i < M; i++)
	  for(int j = 0; j < L; j++)
	    {
	      res.data[i * L + j] = 0;

	      for(int z = 0; z < N; z++)
		res.data[i * N + j] += data[i * N + z] * m.data[z * L + j];
	    }

	return res;
      }

      MTYPE operator +(MTYPE m) const
      {
	MTYPE res(*this);

	for(int i = 0; i < M * N; i++)
	    res.data[i] += m.data[i];

	return res;
      }

      MTYPE& operator +=(MTYPE m)
      {
	for(int i = 0; i < M * N; i++)
	    data[i] += m.data[i];

	return *this;
      }

      MTYPE& operator =(const MTYPE& m)
      {
	for(int i = 0; i < M * N; i++)
	    data[i] = m.data[i];

	return *this;
      }

      MTYPE& operator *=(T val) { for(int i = 0; i < M * N; i++) data[i] *= val; return *this; }

      MTYPE& operator /=(T val) { for(int i = 0; i < M * N; i++) data[i] /= val; return *this; }

      MTYPE operator *(T val) { MTYPE res = *this; for(int i = 0; i < M * N; i++) res.data[i] *= val; return res; }

      MTYPE operator /(T val) { MTYPE res = *this; for(int i = 0; i < M * N; i++) res.data[i] /= val; return res; }

      MTYPE& clear()
      {
	for(int i = 0; i < M * N; i++)
	    data[i] = 0;

	return *this;
      }

      MTYPE& identity()
      {
	if(isSquare())
	  {
	    clear();
	    for(int i = 0; i < M; i++)
	      data[i * M + i] = 1;
	    return *this;
	  }
	else
	  throw PYTHAGORAS_MATRIX_NOT_SQUARE_EXCEPTION;
      }

      constexpr int getRows() { return M; }
      constexpr int getColumns() { return N; }
      constexpr bool isSquare() { return N == M; }

      Matrix<N, M, T> transpose()
      {
	Matrix<N, M, T> res;
	for(int i = 0; i < M; i++)
	  for(int j = 0; j < N; j++)
	    res.data[j * M + i] = data[i * N + j];

	return res;
      }

    protected:
      T data[M * N];
    };
  }
}

#endif //_ICECODE_PYTHAGORAS_MATRIX
